Assessment  of  Translational  Anisotropy  in  Rarefied  Flows 

Using  Kinetic  Approaches  (Preprint) 

D.  C.  Wadsworth  and  N.  E.  Gimelshein*,  S.  F.  Gimelshein'  and  I.  J.  Wysong** 

*ERC,  Inc.,  Edwards  AFB,  CA  93528 
^  University  of  Southern  California ,  Los  Angeles,  CA  90089 
** Propulsion  Directorate,  Edwards  AFB,  CA  93528 

Abstract.  Several  distinct  modeling  approaches  are  tested  for  their  ability  to  resolve  translational  nonequilibrium  phenomena 
that  are  common  to  (but  not  unique  to)  rarefied  gas  dynamic  flowfields.  The  methods  considered  include  DSMC,  BGK,  ES- 
BGK,  Navier- Stokes,  and  an  anisotropic  model  (based  on  a  first  order  expansion  of  the  Boltzmann  equation  for  a  Maxwell  gas 
possessing  an  anisotropic  Maxwellian  distribution  function).  For  the  large  Mach  number  normal  shock  wave,  ES-BGK  shows 
fair-to-good  agreement  with  DSMC,  while  numerical  solution  of  the  anisotropic  model  equations  are  rather  poor.  A  semi- 
analytical  solution  for  the  anisotropic  model  compares  rather  well  with  DSMC  results  for  the  normal  shock  wave,  suggesting 
that  the  equation  set  may  indeed  be  quantitatively  useful  for  modeling  translational  nonequilibrium.  For  the  hypersonic  sphere- 
cone  flowfield,  comparison  of  the  directly  sampled  shear  stress  with  the  Navier-Stokes  (Chapman-Enskog)  and  anisotropic 
value  evaluated  using  the  DSMC  macroscopic  properties  shows  that  the  anisotropic  result  offers  qualitative  improvement  over 
the  NS  value  over  a  broad  range  of  rarefaction. 


INTRODUCTION 

Anisotropy  of  the  molecular  velocity  distribution  function  is  important  for  rarefied  gases  characterized  by  strong 
thermal  non-equilibrium.  Such  a  translational  anisotropy  arises  both  in  compression  flows  dominated  by  shock  waves, 
and  expanding  flows  such  as  nozzle  expansions  and  plumes.  The  need  to  account  for  the  anisotropy  in  nonequilibrium 
gases  was  understood  as  early  as  1867  [1].  Besides  the  transitional  rarefied  flow  regime,  this  anisotropy  may  be 
important  in  turbulent  flow  modeling,  especially  when  the  transition  to  turbulence  needs  to  be  predicted.  The  latter 
one  is  traditionally  approached  with  continuum  methods  based  on  the  solution  of  Navier-Stokes  equations,  which,  as 
was  pointed  out  in  [2],  may  not  be  the  correct  approximate  solution  of  the  Boltzmann  equation  for  this  case.  In  [2],  an 
anisotropic  fluid  seven  equation  set  was  presented  for  the  density,  three  fluid  velocity  components,  and  three  directional 
thermal  kinetic  energies.  This  technique  represents  a  macroscopic  (20  moment)  approach  to  non-equilibrium  flows,  and 
as  such  complements  microscopic,  kinetic  methods,  in  their  consideration  of  flow  nonequilibrium.  A  semi-analytical 
solution  to  this  equation  set  is  also  available  [3]  for  the  normal  shock.  In  the  past,  several  researchers  have  attempted 
to  include  translational  non-equilibrium  in  a  macroscopic  fashion.  Candler  et  al.  [4]  developed  a  multi-temperature 
equation  set  and  compared  numerical  solutions  for  the  normal  shock  wave  problem  with  DSMC  results.  Dogra  et  al.  [5] 
presented  calculations  using  a  multi-temperature  model  for  unsteady  blast- type  problems.  Baganoff  [6]  has  presented 
a  numerical  method  for  Maxwell’s  moment  equations  for  the  normal  shock  problem. 

The  main  objective  of  this  work  is  an  assessment  of  the  accuracy  of  translational  nonequilibrium  predictions 
obtained  with  different  approaches.  In  addition  to  [2]  and  [3],  the  following  numerical  methods  are  used:  the  direct 
simulation  Monte  Carlo  method,  that  is  considered  to  provide  benchmark  solutions;  the  finite  volume  solution  of  model 
kinetic  equations,  and  the  solution  of  Navier-Stokes  equations  as  a  reference  point  in  comparing  macroparameters. 
SMILE  tool  [7]  was  used  for  the  DSMC  computations,  SMOKE  solver  recently  developed  at  ERC  provided  the  model 
kinetic  equation  capability,  and  the  commercial  Fluent  [8]  code  was  applied  to  obtain  Navier-Stokes  results. 

The  DSMC  method  is  considered  the  most  accurate  approach  for  the  present  study,  and  indeed  supplies  the  reference 
data  by  which  the  other  methods  are  evaluated.  A  deterministic  mathematical  equation  set,  however,  can  offer  benefits 
in,  e.g.,  order  of  magnitude  analysis,  which  can  greatly  assist  in  interpreting  the  physical  phenomena  involved.  Two 
model  problems  are  considered  to  assess  nonequilibrium:  the  structure  of  a  normal  shock  wave  over  the  Mach  number 
range  1.2  <  Moo  <  10,  and  the  hypersonic  flow  past  a  blunt  sphere  cone  at  moderate  rarefaction  10-3  <  Kn  <  10-1. 
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NUMERICAL  METHODS 


DSMC 

The  2D/axisymmetric  capability  of  SMILE[7]  was  used  for  the  DSMC  computations.  The  majorant  frequency 
scheme[9]  was  used  to  calculate  intermolecular  interactions.  The  intermolecular  potential  was  assumed  to  be  a 
variable  hard  sphere[10]  with  the  value  of  co  =  0.5  in  the  viscosity-temperature  dependence  in  order  to  simulate  the 
Maxwell  molecule  interaction  used  in  all  other  methods  applied  in  this  work.  For  one-dimensional  normal  shock  wave 
calculations,  specular  surfaces  were  used  to  provide  a  pseudo- ID  setup.  For  2D  calculation  of  an  external  flow,  fully 
diffuse  reflection  was  assumed  for  gas-surface  collisions.  The  total  number  of  simulated  molecules  and  collision  cells 
for  ID  runs  was  approximately  1  million  and  1,000,  respectively,  and  for  2D  runs,  5  million  and  1  million,  respectively. 


ES-BGK 

A  finite  volume  2D/axisymmetric  code  SMOKE  developed  at  ERC  has  been  used  to  solve  the  BGK  and  ES-BGK 
equations.  SMOKE  is  a  parallel  code  based  on  conservative  numerical  schemes  developed  by  L.  Mieussens  [11].  A 
second  order  spatial  discretization  is  used  along  with  implicit  time  integration.  1,000  spatial  cells  were  used,  with 
specular  boundaries  set  to  model  ID  flow  in  the  longitudinal  direction.  The  convergence  on  the  velocity  grid  was 
studied,  with  the  number  of  (v,y,z)  velocity  points  ranging  from  (30,20,20)  to  (60,40,40).  The  results  for  the  latter  one 
are  presented  below. 


Navier-Stokes 

Calculations  are  performed  with  a  standard  commercial  CFD  code,  Fluent  [8].  A  density-based  Navier-Stokes 
solution  was  obtained  for  a  spatially  uniform  grid  with  1,000  cells;  the  Maxwell  molecule  power-law  viscosity- 
temperature  dependence  was  used,  fl  T  in  order  to  match  the  viscosity  used  in  the  kinetic  methods. 


Anisotropic  Method 

Kliegel  [2]  presents  the  continuity,  momenta,  and  directional  kinetic  energy  equations  derived  from  a  first  order 
expansion  of  the  Boltzmann  equation  for  a  Maxwell  gas  possessing  the  anisotropic  or  multi-temperature  near¬ 
equilibrium  velocity  distribution  function, 

U(C)  =  (2*K)-^(W,)-'/;e*p  (-  +  mj)  '  (1) 

with  the  mean  temperature  T  =  (7\  +  T2  +  T^)/3.  The  equation  set  is  similar  to  that  given  in  [4]. 


Solution  Method 

The  semi-analytic  normal  shock  wave  solution  is  available  [3]  in  tabular  form  for  a  range  of  Mach  numbers.  For  the 
numerical  solution,  the  public  domain  FiPy  [12]  finite  volume  partial  differential  equation  solver  is  used  to  integrate 
equations  (31)-(33)  of  [2]. 


RESULTS 

Normal  shock  structure 

Figure  1  compares  the  calculated  profiles  of  normalized  temperature  and  velocity  for  Mach  1.2.  For  this  relatively 
weak  shock  wave,  all  results  are  comparable,  as  expected.  Numerical  and  semi-analytic  solutions  to  the  anisotropic 
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equations  are  found  to  be  nearly  identical,  and  only  the  numerical  solution  is  presented  in  this  figure.  Note  that  not  only 
macroparameters  obtained  by  different  approaches  are  close  for  this  case,  but  also  the  velocity  distribution  functions. 

Figures  2  and  3  compare  the  profiles  for  Mach  10.  Under  these  conditions,  the  numerical  solution  of  the  equations 
of  Kliegel  [2]  (not  shown)  was  found  to  be  poor.  Candler  et  ah  have  presented  results  for  Mach  1 1 ,  finding  similar 
issues,  and  suggest  that  retaining  additional  terms  in  the  expansion  could  improve  the  robustness  of  the  equation  set. 
The  semi-analytic  solution  of  the  equations  [2]  is  in  very  good  agreement  with  the  DSMC  solutions.  Note  also  that  the 
semi-analytic  solution  does  predict  a  small,  about  1%,  overshoot  in  the  mean  temperature,  that  is  typically  observed  in 
the  DSMC  method  (see,  for  example,  [13]).  This  overshoot,  although  somewhat  larger  than  in  DSMC,  is  also  predicted 
in  the  solution  of  the  ES-BGK  equation.  The  overshoot  is  not  seen  in  the  solution  of  the  BGK  and  the  Navier- Stokes 
equations.  As  expected,  the  latter  ones  are  characterized  by  much  thinner  shock  than  the  other  approaches. 

It  is  useful  to  consider  equation  (76)  of  [2]  for  the  longitudinal  temperature  T\ ,  obtained  by  combining  the  continuity 
and  longitudinal  momentum  equations, 


Too 


U\ 

Uoo 


yMl  + 1 


(2) 


where  y  is  the  gas  ratio  of  specific  heats,  and  the  subscript  °o  indicates  upstream  values.  Such  an  analytic  expression 
for  T\  as  a  function  of  longitudinal  flow  velocity  u\  serves  as  an  additional  indicator  of  the  numerical  accuracy  of 
obtained  solutions. 


TABLE  1.  Maximum  value  of  normalized  T\  for  different  approaches. 


Exact 

DSMC 

Anisotropic  Semi-analytic 

Anisotropic  Numerical 

ES-BSK 

BGK 

1.3227 

1.3206 

1.3226 

1.2579 

1.3227 

1.3227 

This  equation  predicts  overshoot  in  T\  arises  for  Moo  >  y/9/5  ~  1.34.  A  smaller  (few  percent)  overshoot  in  the 
mean  temperature  also  arises  (see  Fig.  3).  As  seen  in  Table  1,  where  the  value  of  ^  is  shown  at  its  maximum 
at  ui/u-oo  =  0.5,  the  semi-analytical  solution  recovers  the  overshoot  in  temperature  seen  in  the  DSMC  and  ES-BGK 
results.  Note  that  the  DSMC  value  is  slightly  lower  than  the  exact  value  given  by  equation  2  due  to  statistical  scatter  and 
a  finite  time  step.  The  numerical  solution  of  the  anisotropic  equations  was  found  to  underpredict  the  peak  T\  values,  in 
agreement  with  the  results  of  [4];  the  longitudinal  temperature  profile  obtained  from  (2),  however,  was  found  to  yield 
an  accurate  profile.  This  suggests  that  the  modifications  to  at  least  the  longitudinal  energy  transport  equation,  along 
with  more  sophisticated  numerical  integration  schemes,  may  be  necessary  to  capture  the  strong  nonequilibrium  in  this 
problem. 


FIGURE  1.  Mach  1.2  normal  shock 
macroparameters  for  different  numeri¬ 
cal  methods. 


FIGURE  2.  Mach  10  normal  shock 
flow  velocities  for  different  methods. 


FIGURE  3.  Mach  10  normal  shock 
mean  temperatures  for  different  meth¬ 
ods. 


Figure  4  compares  the  directly  sampled  velocity  distribution  functions  obtained  from  DSMC,  ES-BGK,  and  BGK 
solutions  at  the  position  of  maximum  T\  overshoot  for  the  Mach  10  shock.  The  dotted  lines  show  the  corresponding 
equilibrium  distribution  evaluated  using  local  macroscopic  properties.  Although  all  major  macroparameters  obtained 
by  different  approaches  are  very  similar  for  a  given  value  of  longitudinal  temperature,  the  distribution  functions  are 
noticeably  different,  especially  for  the  BGK  solution.  The  dissipation  of  the  peak  corresponding  to  the  free  stream  flow 
is  almost  complete  in  the  region  near  the  maximum  7j ,  while  for  the  other  two  approaches  the  distribution  function  is 
clearly  bimodal.  Note  that  similar  bimodal  behavior  was  observed  for  BGK  upstream  of  the  maximum  7j. 
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The  deviation  of  the  model  distributions  from  the  directly  sampled  value  can  be  evaluated  at  each  point  through 
the  shock  and  used  as  a  measure  of  the  relative  improvement  they  may  provide.  Here,  the  deviation  is  defined  as  a 
simple  density-weighted  value,  S  =  f  \ fosMC  —  f\dc\.  Figure  5  shows  the  calculated  profiles.  In  this  figure,  all  three 
distribution  functions  are  based  on  the  macroparameters  obtained  from  the  Mach  10  DSMC  computation,  so  they  can 
be  compared  directly.  For  this  admittedly  arbitrary  parameter,  no  model  gives  a  clearly  superior  approximation  to 
the  actual  DSMC  distribution  function.  Generally,  detailed  comparison  of  distribution  functions  obtained  by  different 
approaches  showed  that  the  solution  of  the  ES-BGK  equation  provides  the  best  approximation  to  the  distribution 
functions  obtained  in  DSMC  in  terms  of  the  shape  and  capturing  the  bimodal  behavior.  Note  also  the  use  of  the  DSMC 
distribution  functions  as  the  reference  point  is  complicated  by  the  fact  that  they  were  found  to  depend  on  the  scattering 
law  [14]. 


FIGURE  4.  Longitudinal  velocity  distribution  function  di¬ 
rectly  sampled  in  DSMC,  ES-BGK,  and  BGK,  at  the  point 
where  T\  is  maximum. 


FIGURE  5.  Deviation  of  Maxwellian,  Chapman-Enskog, 
and  anisotropic  longitudinal  velocity  distribution  functions 
from  DSMC. 


Moments  of  the  distribution  function  are  of  more  direct  interest,  with  higher  order  moments  providing  a  critical 
test  of  the  relative  accuracy  of  the  models.  Figure  6  shows  profiles  of  the  4th  order  moment  (c\)  inside  the  Mach 
10  shock,  obtained  by  four  different  approaches.  The  semi-analytic  anisotropic  model  agrees  relatively  well  with  the 
directly  sampled  DSMC  results,  generally  providing  noticeably  better  agreement  than  the  solution  of  the  ES-BGK 
equation.  Interestingly,  the  maximum  value  of  the  4th  order  moment  obtained  by  the  continuum  approach  (NS)  is 
close  to  DSMC. 


FIGURE  6.  Profiles  of  4th  velocity  moment  ratio  obtained  FIGURE  7.  Profiles  of  the  ratio  of  two  3rd  order  moments 

with  four  different  numerical  methods.  obtained  with  four  different  numerical  methods. 
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An  important  indicator  of  the  ability  of  a  numerical  approach  to  capture  flow  nonequilibrium  is  the  ratio  of  two  3rd 
order  moments,  {c\) / {c\c2).  The  numerator  is  the  normal  energy  flux  and  the  denominator  is  the  normal  total  energy 
flux.  This  ratio,  computed  with  four  different  approaches,  is  plotted  in  Fig.  7.  Note  that  the  near-equilibrium  parts, 
both  downstream  and  upstream,  where  the  ratio  becomes  numerically  unstable,  are  not  shown  here.  The  solution  of 
the  Navier-Stokes  equations  provides  a  constant  value  of  0.6  for  this  ratio,  which  does  not  depend  on  the  location 
inside  the  shock  (0.6  is  in  fact  the  equilibrium  value).  In  the  DSMC  solution,  the  ratio  varies  from  about  0.75  in  the 
upstream  region  to  about  0.38  in  the  downstream  region.  The  semi-analytic  anisotropic  model  agrees  well  with  the 
DSMC  results,  whereas  the  ES-BGK  solution,  being  qualitatively  similar,  noticeably  overpredicts  them. 


Sphere  cone  flow 


The  flowfield  surrounding  a  vehicle  in  hypersonic  high  altitude  flight  is  a  more  realistic  and  interesting  test  of  the 
applicability  and  accuracy  of  any  nonequilibrium  model.  Here,  the  vehicle  geometry  is  a  10  degree  half  angle  cone 
with  a  spherical  nosetip  of  radius  rn,  and  overall  length  /  =  10r„.  The  vehicle  speed  is  4.5 km/s  (Mach  number  on  the 
order  of  15-20),  and  the  freestream  density  is  varied  to  obtain  Kni  ~  10-3, 10-2, 10-1.  The  property  examined  here  is 
the  shear  stress  evaluated  by  different  approaches. 

Figure  8  shows  contours  of  the  shear  stress  for  Kn  =  10-3.  The  DSMC  values  are  directly  sampled  {c\C2).  The 
anisotropic  values  are  evaluated  using  the  DSMC  macroscopic  values  and  equation  (66)  of  [2], 


du\  duo 

^2  "5  ~  H-  T\  -r— - 
0X2  dx\ 


The  Navier-Stokes  values  are  evaluated  using 


Tl2  =  M 


du\ 

dx2 


du2 

dx\ 


(3) 

(4) 


The  anisotropic  model  reduces  to  this  form  for  T  =  7}  and  using  the  definition  of  the  viscosity  for  the  Maxwell  model, 

For  these  near-continuum  conditions,  all  results  are  comparable,  as  expected.  The  difference  between  shear  stress 
values  obtained  with  different  approaches  is  within  a  few  percent  in  the  bow  shock,  and  only  slightly  larger  in  the 
boundary  layer.  The  biggest  difference,  up  to  a  factor  of  two,  is  observed  in  the  expansion  region  immediately  after 
the  cone  rear  edge.  The  DSMC  values  are  lower  in  this  region,  and  those  obtained  for  the  anisotropic  model  and  NS 
are  close. 


x,  m 


Stress 


y 


60.00 

56.47 

52.94 

49.41 

45.88 

42.35 

38.82 

35.29 

31.76 

28.24 

24.71 

21.18 

17.65 

14.12 

10.59 

7.06 

3.53 

1.87 

0.91 


2 


X,  m 


Stress 

13.000 

10.062 

7.788 

6.028 

4.666 

3.611 

2.795 

2.163 

1.674 

1.296 

1.003 

0.776 

0.601 

0.465 

0.360 

0.279 

0.216 

0.167 

0.129 

0.100 


FIGURE  8.  Shear  stress  contours,  Kn  =  FIGURE  9.  Shear  stress  contours,  FIGURE  10.  Shear  stress 

10-3:  DSMC,  Anisotropic,  NS  Kn  =  10-2:  DSMC,  Anisotropic,  NS  contours,  Kn  =  10-1:  DSMC, 

Anisotropic,  NS 


Figure  9  shows  contours  of  the  shear  stress  for  Kn  =  10-2.  For  this  transitional  rarefaction  case,  the  anisotropic 
results  show  qualitative  improvement  relative  to  the  NS  values  throughout  the  shock  layer.  Figure  10  shows  contours 
of  the  shear  stress  for  Kn  =  10-1.  For  this  rarefied  case,  the  anisotropic  results  show  even  a  better  improvement  than 
for  Kn  =  10-2  (factor  of  3-5),  especially  in  the  shock  region  and  in  the  boundary  layer. 
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CONCLUSIONS 


Several  kinetic  (the  DSMC  method,  the  solution  of  BGK  and  ES-BGK  equation)  and  continuum  (the  solution  of 
Navier-Stokes  equations  and  the  anisotropic  method  of  Kliegel)  approaches  are  tested  in  this  work  for  their  ability 
to  accurately  predict  translational  nonequilibrium  effects.  The  DSMC  solutions  are  considered  benchmark,  and  the 
other  approaches  are  compared  with  them  on  a  microscopic  level  (velocity  distribution  functions)  and  macroscopic 
level  (major  macroparameters  and  higher  moments).  The  problems  under  consideration  are  the  structure  of  the  normal 
shock  wave,  both  weak  and  strong,  and  a  hypersonic  flow  over  a  blunted  cone. 

All  used  numerical  approaches  are  found  to  perform  well,  both  macro  and  microscopically,  for  a  Mach  1 .2  shock 
wave,  where  the  deviation  from  equilibrium  is  relatively  small.  For  a  Mach  10  shock,  significant  differences  from 
DSMC  were  observed  for  both  kinetic  and  continuum  methods.  While  the  solution  of  the  ES-BGK  equation  provides 
better  qualitative  agreement  with  DSMC  in  terms  of  the  longitudinal  velocity  distribution  functions,  the  anisotropic 
model  was  found  to  perform  better  for  both  lower  and  higher  order  moments  of  the  distribution  function. 

In  a  hypersonic  flow  over  a  blunted  cone,  the  expression  for  the  shear  stress  used  in  the  anisotropic  model  are 
shown  to  offer  significant  improvements  as  compared  to  the  conventional  Navier-Stokes  expression.  Generally,  the 
anisotropic  equation  set  may  be  useful  for  flows  with  moderate  flow  nonequilibrium,  and  may  have  some  potential 
in  addressing  problems  associated  with  the  transition  to  turbulence.  For  severely  nonequilibrium  flows,  higher  order 
terms  may  have  to  be  included  in  the  anisotropic  equations.  The  relative  importance  of  the  assumption  of  Maxwell 
molecules  used  in  the  anisotropic  equations  is  still  to  be  examined. 
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